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We study the evolution of inhomogeneous spherical perturbations in the universe 
in a way that generalizes the spherical top hat collapse in a straightforward manner. 
For that purpose we will derive a dynamical equation for the evolution of the density 
\ contrast in the context of a Lemaitre-Tolman-Bondi metric and construct solutions 
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with and without a cosmological constant for the evolution of a spherical perturbation 



with a given initial radial profile. 
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The formation of large scale structure in the universe is one of the most promising probes 
00 ■ 

that will be used to determine fundamental cosmological parameters in future extra-galactic 
surveys, such as the imminent Dark Energy Survey (DES) [1], the Large Synoptic Survey 
Telescope (LSST) [2] and the Euclid survey j^]. 

The basic picture for large scale structure formation is that small perturbations gener- 



ated by quantum fluctuations in the inflationary epoch grow in the dark matter to form 
gravitational potential wells where the baryons later fall into. The linear part of thispro- 



Q. 



cess for dark matter is well understood and described by linear perturbation theory [4 
However, the nonlinear stages are not amenable to these perturbative methods even in the 
case when only dark matter is present. One then usually resorts to large numerical N-body 
simulations in order to obtain, e.g. the mass distribution of large dark matter haloes. Un- 
fortunately, these simulations are very costly and it is desirable that a simple, approximate 
semi-analytical approach could be used to estimate the properties of dark matter haloes. 



One such model is the so-called spherical collapse model [6]. In this extremely simple 
model, a spherically symmetric region of homogeneous overdensity, called "top-hat" density 
profile, evolves inside a homogeneous expanding Universe. Symmetry arguments show that 
one can regard the overdense region as a mini-universe of positive curvature. Hence the 
Raychaudhury and continuity equations can be used to evolve the density and radius of the 
spherical region j3> 8|. When combined with the Press-Schechter theory [9] this framework 
provides a statistical basis for structure formation from which the number density of dark 
matter haloes can be estimated. 

The spherical collapse model has been recently used to study structure formation in 
models with a simple Yukawa-type modification of gravity jlO|, in the so-called f(R) models 
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11 1 . in braneworld cosmologies 121 ]. in models which allow for dark energy fluctuations 



^23] and in the so-called chameleon models |24j |. The purpose of this work is to relax 
the assumption of a top-hat profile and study the nonlinear evolution of an inhomogeneous 
spherical perturbation in a fully relativistic framework. 



II. EVOLUTION OF AN INHOMOGENEOUS SPHERICAL PERTURBATION 

In the case of an inhomogeneous perturbation one can not use a Friedmann-Robertson- 
Walker (FRW) metric which is valid only for a homogeneous matter distribution. The most 
general spherically symmetric metric is given by: 

ds 2 = e A ^dt 2 - e B ^dr 2 - R 2 (r, t)dtf (1) 

where t is the cosmic time, r is the comoving radial coordinate and dQ is the solid angle 
element. The metric is determined by three functions: A(r, t), B(r, t) and R(r, t), where this 
last function is known as the areal radius since the area of a surface with a given time and 
comoving radial coordinates is given by S = AuRirA) 2 . For a nice discussion on spherically 

s ymm e tr ic i „ho m o g eneous m ode 1 ssee,e.,,chap t e rl Sina. 

We assume a single perfect fluid with energy-momentum tensor given by 

T a b = (p + p)u a u b -p8\ (2) 

where p(r, t) and p(r, t) are the energy density and pressure of the fluid and u a is its 
4— velocity, with u a u a = 1. In a comoving reference frame u = (e~ A ^ 2 , 0). 
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In this metric, Einstein equations with a cosmological constant, G a f, — Ag a b = 8nGT a i n 
can be written as: 



• / R 2 BR \ d ( R" R' 2 B'R'\ 1 „ . . 



, i R R 2 AR \ r> / i?'^ v4'i?'\ 1 



^ [ 4^ - 2^ + 2^ + 2fi + - Si? 1 - 
4 I R R R 

e~ B ( R" A'R' B'R' 



(R" A' R' R' R' \ 
4^- + 2^- - 2^- + 2A" + A' 2 - A'B' = -8nGp + A (6) 
R R R J 

where dots and primes refer to partial derivatives with respect to time and space, respectively. 

Conservation of the energy-momentum tensor, V a T a b = 0, in this spacetime results in the 

following equations: 

fi + 4^ + 2^ = (7) 
R p+p K 1 

A' + 2^— = Q. (8) 
P + P 

Equations and combined with Eq. (J3J), can be written as two conservation equa- 
tions for the so-called active mass m(r, t): 

2m'(r, t) 



R?R' 

2m(r, t) 

R 2 R 

with 



87rGp (9) 
-8nGp (10) 



2m(r, t) = e~ A RR 2 - e~ B RR' 2 + R- ^AR 3 . (11) 
We will focus on pressureless dark matter, in which case p = and 

m(r,t) = =>- m = m(r). (12) 



In this case the metric takes the form of the Lemaitre-Tolman-Bondi (LTB) metric 26M28]: 



ds 2 = dt 2 Krrdr 2 - R 2 (r, t)dn 2 (13) 

1 + f(r) 
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where the curvature f(r) is determined by m(r) and R 

m = &-*?p.- Ur*. (14) 

Models with large inhomogeneities described by a LTB metric, such as the so-called "Hub- 
ble bubble" model, have been used recently as an alternative explanation to the apparent 



acceleration of the universe (see, e.g 29l432j|). Our goal here is to study the cosmological 
evolution of such an inhomogeneity. 

The dynamical equations for the LTB model can be written as a generalization of the 
Friedmann equation 

R 2 _ 2m(r) f(r) A 

and a generalization of the "acceleration" equation 

R m(r) 



(16) 



R R 3 ' 

which is independent of the curvature, just like the FRW case. These equations reduce to 
the usual FRW equations if we set 

R(r,t) = a(t)r (17) 

and 

f(r) = -kr 2 , (18) 

where k is the spatial curvature. 

To describe the evolution of the inhomogeneous density perturbation we define the density 
contrast 

s^eMzm (19) 

where p(r, t) is the energy density inside the perturbation and p(t) is the energy density of 
the unperturbed expanding background. 

The dynamical equation for 5(r, t) can be deduced from Eq. ([7j), which for a pressureless 
fluid becomes 

p + 3hp = (20) 

where h = b/b and b = (R 2 R') l l 2 , together with the expanding background continuity 
equation: 

p- + 3Hp = (21) 
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where H = a/a. Taking the derivative of Eq. (TL9]) with respect to time we obtain 

P = p(l + 5) + p5, (22) 

which can be rewritten as 

5 = 3 (H - h) (1 + 5) (23) 



where we have used Eqs. (120]) and fl2T|) . 

Deriving Eq. (123]) with respect to time we obtain 



5 = 3(H-h) (l + 5)+3(H-h)5. (24) 



If we now combine Eqs. ([3])-(j6]) with the equations 



2 2 

and 



H = -Ih 2 + ^ (25) 



h = 4nGp + ^- ( ^ + 4 ) + A - 3h ' 2 ( 26 ) 



3R \R R 

Eq. (124"]) becomes the nonlinear differential equation for the evolution of the density pertur- 
bation: 



4 5 2 A „ „ s 2 r 1 



5 + 2#5-47rGp5(l + 5)--— = _(i + J 5) + _ 



9, i? 



1 + S) . (27) 



Note that for a homogeneous density perturbation with a top-hat profile the areal radius 
is R = ra. In this case, the last term on the right hand side of Eq. (1271) vanishes, reducing 



this equation to the well known evolution equation (see, e.g. 19]) in the case of a single 
pressureless fluid. 

Eq. (T2T]) generalizes the evolution of 5 for spherical perturbations with arbitrary initial 
radial profiles. Note, however, that for the simple case of pressureless dark matter, the 
solution of this equation can be obtained directly from Eqs. (120]) and (9), which is what we 
will do in the following. 

Consider initially a radial density profile expanding with the background 

Ri = rai (28) 
which is characterized by a time-independent function g{r) 

p l (r)=p l {l + 5 l {r)} = p l g{r) . (29) 
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Since Eq. (9) in LTB is time-independent, we can set 

AnGR 2 R' iPi (r) = A7rGR 2 R'p . (30) 



Furthermore, since p — P fl3 an d P — P {1 + <5} we obtain 

R 2 R 



Hr,t)=a(tf^£-1. (31) 



Defining a dimensionless areal radius 



R = -, (32) 

r 



Eq. (}3~Tj) becomes 

g( r ,t) = a(t) 3 - g,^ (33) 
R 2 R + rR' 

which is, of course, the formal solution to the nonlinear evolution equation Eq. (12T1) . 

Before we find the evolution of the density contrast we need to solve the dynamical 
equations of the background FRW metric for a(t) and of the LTB metric for R(r,t). 



III. SOLUTIONS OF THE LTB MODEL FOR A GIVEN INITIAL DENSITY 

PROFILE 

We consider initially solutions of the LTB model applied to an inhomogeneous density 
perturbation with an initial density profile embedded in a Einstein-de Sitter and ACDM 
background. 



A. Einstein-de Sitter background 

Our first goal is to study the evolution of an initial spherically symmetric perturbation 
in a dark matter dominated universe with an Einstein- de Sitter (EdS) background in the 
context of General Relativity. 

We assume that at an arbitrary initial time t = ti, the whole universe follows the back- 
ground expansion with R(r,ti) = a(ti)r, but has an initial small density perturbation spec- 
ified by a profile g(r): 

p(r,U) = p(ti)g(r) (34) 
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where p{tj) is the initial homogeneous background density. 
Since m(r) is time-independent it can be computed as: 

m(r) = AixGalpi I dr'r' 2 g(r'). (35) 
Jo 

Defining 

h(r) = 4 f dr'r' 2 g{r') (36) 
r Jo 

we can write 

m ( r ) = —^-p~or 3 h(r), (37) 

where the background density today is p Q = afpi. 
Using time in units of the Hubble age today 

t = tH (38) 

and distance in units of the comoving radius 

R = - (39) 
r 

we will solve the acceleration equation 

^ = -\^R' 2 Kr) (40) 
for a given value of r and initial conditions 



= + (1 - fio) (41) 

and Ri = a«. 

Hence, given the initial conditions, the evolution of the inhomogeneous perturbation is 
completely determined by its initial profile. As an example of our method, we choose a 
Gaussian profile with width given by r c with a sharp cutoff at r = r c characterizing the 
comoving size of the perturbation: 

g{r) = l + 5 l e- r2 l r "Q{r c -r) (42) 

where Q(x) is the usual Heaviside function. We have started our evolution from = 1CT 5 
and chose 5i = 10" 4 552 , which results in the collapse today (a = 1 when t = 0.667) of the 
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center of the perturbation. Exact solutions for the dynamical LTB equations in this case 
exist in a parametric form [25] and we used them to check our numerical method. 

In Fig.(CQ) we show the time evolution of R(r, t) for different values of the comoving radius 
r in units of r c . We can see that different shells collapse at different times, with the outer 
shells collapsing later. Hence, there is no shell crossing in this expected since the 

density profile is decreasing radially. 




0.8 



FIG. 1: Areal radius R(r,t) as a function of time t for r/r c = (solid line), 0.2 (dashed line) and 0.5 
(dot-dashed line). 



Interestingly, the presence of a density perturbation in our example of an EdS universe, 
even when initially localized, makes the whole universe to eventually collapse. For instance, 
a region with a comoving radius of r = 2r c collapses at t — 34.1. Therefore, the expansion 
rate of the universe at this radius is slightly different from the background EdS. 

In Fig.(T5]) we show the initial instants of the evolution of the profile of the density contrast 
5(r, t). The behavior of the perturbation near the sharp boundary is a consequence of 
numerical instabilities that do not appear in the exact solution. 

When one is still in the linear regime, the shape of the profile does not change significantly. 
For times close to the collapse of the center, the density profiles approach Gaussian profiles 
with different widths, as shown in Fig. ([3]). 
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FIG. 2: Density contrast 5(r, t) as a function of comoving radius in units of r c for the initial contrast 
(solid line), density contrast at i = 10~ 7 (dashed line) and t = 10~ 6 (dot-dashed line). 
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FIG. 3: Density contrast 5{r, t) as a function of comoving radius in units of r c for t = 0.5 (dashed line), 
density contrast at t = 0.60 (dot-dashed line) and t = 0.61 (solid line). 

B. ACDM background 



We now proceed to study the evolution of an inhomogeneous spherical density perturba- 
tion in a background with dark matter and cosmological constant, where ~p{t) = ^ + Pa 
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and n° m + n° A = 1. 

As the cosmological constant is not perturbed, it only affects the behavior of the density 
contrast through changing the evolution of a (i) and R(r,t). We solve the LTB equations 
with a cosmological constant, assuming that at an arbitrary initial time ti = 0, the whole 
universe follows the background expansion with R(r,ti) = a(ti)r but it has an initial small 
density perturbation specified by the profile g(r). 

The function mass m(r) now becomes 

m (r) = l -Hlnlr*h (r) . (43) 

Using time in units of the Hubble age today and distance in units of the comoving radius, 
we will solve the acceleration equation 

d 2 R 1 



dt 2 2 

for a given value of r and initial conditions 



-ST m R- 2 h (r) + Q A R (44) 




3* + nw (45) 

and Ri = a$. In our example below we will use = 0.23 and Q° A = 0.77, in which case 
t = 1.037 today. 

We should point out that also in this case there are exact solutions in terms of implicit 
Weierstrass-p functions for which the implementation has to be carried semi-analytically 



33|. Since this implementation is perhaps as intricate as solving the equations numerically, 
we have opted for the latter procedure, for simplicity. 

In Fig. (j3J) we show that the time evolution of R(r, t) has the same qualitative behavior as 
in the EdS case, except that from a certain radius on (around r = 0.7r c in our example) the 
effect of the cosmological constant halts the collapse. This is in agreement with the recent 



analysis of 34], where it was found a dividing shell separating expanding and collapsing 
regions in general models. 

The initial evolution of perturbations in matter is shown in Fig. (IS]). Again we show the 
initial instants of the evolution of the profile of the density contrast 5(r, t). The small bump 
around r = r c is due to numerical instabilities and appears because of the sharp boundary 
in the density profile. Notice that the perturbation profile rapidly evolves to a Gaussian in 
the nonlinear regime, as can be seen in Fig. ([6]). 
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FIG. 4: Areal radius R(r,i) as a function of time t for r/r c = 0.01 (solid line), 0.5 (dashed line), 0.7 
(dot-dashed line) and 0.75 (dotted line). 
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FIG. 5: Density contrast 5(r, i) as a function of comoving radius in units of r c for the initial contrast 
(solid line), density contrast at t = 10 -7 (dashed line) and t = 10~ 6 (dot-dashed line). 

IV. CONCLUSION 



In this work we have implemented a simple generalization of the spherical top-hat collapse 
model including a general initial profile for the perturbation. We compute the nonlinear 
evolution of the density perturbation. We chose as an example a simple Gaussian profile 
with a sharp cutoff that separates the perturbed region from the background. We showed 
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FIG. 6: Density contrast 5(r,i) as a function of comoving radius in units of r c at t = 0.5 (solid line), 
t = 0.6 (dashed line) and t = 0.7 (dot-dashed line). 



that the density perturbations evolve differently in an EdS and a ACDM background, where 
in the latter case there is a dividing shell between an expansion and a contraction region 
inside the perturbation. One should notice that from the usual definition of the quantity 
S c (z), which is the linearly evolved density contrast with an initial condition such that the 
nonlinear collapse occurs at a redshift z, that its value in principle depends only on the 
background cosmology, being independent of the density profile. 

We plan to extend our method to more realistic models including interacting fluids with 
pressure, where no closed form solutions can be obtained. 
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